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Lattice QCD thermodynamics with Wilson quarks 

Shinji Ejiri*) 

Physics Department, Brookhaven National Laboratory, 
Upton, New York 11973, USA 

We review studies of QCD thermodynamics by lattice QCD simulations with dynamical 
Wilson quarks. After explaining the basic properties of QCD with Wilson quarks at finite 
temperature including the phase structure and the scaling properties around the chiral phase 



energies. 



transition, we discuss the critical temperature, the equation of state and heavy-quark free 

< 

81. Introduction 

<N . 

In order to remove theoretical uncertainties in the analysis of data from heavy- 
ion collision experiments, first principle calculations by the lattice QCD are indis- 
pensable and various interesting results have been already reported. However, most 
lattice QCD studies at finite temperature (T) and chemical potential (/x g ) have been 
4h ■ performed using staggered-type quark actions with the fourth-root trick of the quark 

determinant so far. Therefore, studies by the other actions such as Wilson-type quark 
actions are necessary to estimate systematic errors due to lattice discretization. 

In this report, we want to highlight the QCD thermodynamics by numerical 
simulations with dynamical Wilson quarks. A systematic study of the QCD ther- 
r^- ■ modynamics has been done by the CP-PACS Collaboration using the Iwasaki (RG) 

improved gauge action and the 2 flavor clover improved Wilson quark action sev 



eral years ago. 1 )' 2 ) Recently the WHOT-QCD Collaboration restarted the study by 
the same action adopting new technical developments. 3 ^ 5 ) In this review, we first 
explain in Sec. [2] the phase structure of 2 flavor QCD with Wilson quarks in the 
simulation parameter plane (j3,K). Then, the universality class of the chiral phase 
transition is discussed in Sec. [3l An estimation of the critical temperature in the 
chiral limit is given in Sec. 01 The equation of state (EoS) is discussed in Sec.EJ high- 
lighting fluctuations at finite density, related to the physics of the possible critical 
point in the (T,n q ) plane., and also results of heavy-quark free energies are shown 
in Sec. [H 

§2. Phase structure of QCD with Wilson quarks 

The lattice QCD with Wilson-type quarks is known to have a complicated phase 
structure due to the explicit violation of chiral symmetry and due to the existence 
of the parity-flavor broken phase (Aoki phase). 6 ) Therefore, a systematic study 
surveying a wide range of the parameter space is required to determine appropriate 
simulation parameters. Moreover, lattice artifacts are large on coarse lattices used in 
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most finite temperature simulations when the standard plaquette gauge action and 
the standard Wilson quark action are used, For example, unexpected strong phase 
transition is observed at intermediate quark masses in 2 flavor QCD. Therefore, we 
have to improve the lattice action to reduce these lattice artifacts. 

Figure [1] (left) is the phase diagram obtained by the CP-PACS Collaboration for 
2 flavor QCD with the RG-improved gauge action combined with the clover-improved 
Wilson quark action. ^ The 2 flavor lattice QCD has two simulation parameters K 
and (5. K is the hopping parameter in the quark action and (3 is Q/g 2 in the gauge 
action. The solid line K C (T = 0) is the location of the chilal limit. The pion 
mass decreases as K increases from small K, and vanishes on the line K C (T = 0) 
at zero temperature. In the region above K C (T = 0), a parity-flavor symmetry of 
the Wilson-type quark action is broken spontaneously. 6 ) At zero temperature, the 
boundary of the parity-flavor broken phase is known to form a sharp cusp touching 
the free massless fermion point K = 1/8 at /3 = oo. However, the region above 
K C (T = 0) is an unphysical parameter space, and we usually perform simulations 
below K c . The region below K C (T = 0) corresponds to the physical QCD with 
1/K — 1/K C being proportional to the quark mass m q . 

At finite temperature, the parity-flavor broken phase retracts from the large (5 
region. The colored region with the boundary K C (T > 0) in Fig. [1] (left) is the 
parity-flavor broken phase at finite temperature for the temporal lattice size Nt = 4. 
When Nt is fixed, the temperature T = (A^a) -1 becomes higher as j3 is increased, 
since the lattice spacing a becomes smaller. The dashed line Kt is the pseudo-critical 
line separating the hot and cold phases for Nt = 4. The region to the right of Kt 
(larger (3) is the high temperature quark-gluon plasma (QGP) phase, and that to 
the left (smaller (5) is the low temperature hadron phase. The crossing point of the 
K C (T = 0) and the Kt is the chiral phase transition point. 

As shown in this figure, the line that the pion mass vanishes at T > (K C (T > 
0)) runs along the line of the chiral limit in the low temperature phase, while the 
K C (T > 0) line bends sharply at the chiral phase transition point and goes to the 
unphysical region above K C (T = 0). This is consistent with the picture that the 
massless pion, i.e. the Goldstone boson associated with spontaneous chiral symmetry 
breaking, appears only in the cold phase. 

§3. 0(4) scaling relation 

Next, we discuss the natures of the chiral phase transition in the critical region. 
The order of the phase transition is expected to be second order for 2 flavor QCD and 
first order for 3 flavor QCD in the chiral limit. The confirmation of this expectation 
is an important step toward the clarification of the QCD transition in the real world. 
In this section, we discuss the scaling behavior of 2 flavor QCD. 

When the chiral transition of 2 flavor QCD is second order, the transition is 
expected to be in the same universal class with a 3-dimensional 0(4) spin model. 
With the identifications M ~ for the magnetization, h ~ m q a for the external 

magnetic field and t ~ (3 — (3 c t for the reduced temperature, where /3 c t is the chiral 
transition point, we expect the same scaling behavior as the 0(4) spin model. 
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Fig. 1. Left: Phase diagram for the RG-improved gauge and clover-improved Wilson quark actions 
on an Nt = 4 lattice. Right: 0(4) scaling relation of the chiral condensate. 1 ^ 



The 0(4) scaling was first tested with staggered fermion by Karsch and Laer- 
mann. 7 ) The study was extended to a wider range of the quark mass and lattice sizes 
by several groups. 8 ^ 11 ) However, an agreement of the critical exponents between the 
0(4) spin model and QCD with 2 flavors of staggered quarks has not been obtained. 

For the case of Wilson quarks, Iwasaki et al. 12 ) investigated the scaling relation, 



M/h 1 / 5 = f{t/h l " M ) 



1/PS\ 



(3-1) 



using the standard Wilson quark action coupled to RG-improved gluon. They identi- 
fied the subtracted chiral condensate defined by an axial Ward-Takahashi identity, 13 ) 
(\P\P)sub = 2m q aZ ^2 x (tt(x)tt(0)), as the magnetization of the spin model. Here, the 
quark mass m q is defined by an axial vector Ward-Takahashi identity, 13 ) m^ WI , and 
the tree-level renormalization coefficient Z = (2K) 2 was adopted. They found that 
the scaling relation Eq. (|3-ip is well satisfied with the critical exponents and the 
scaling function of the 0(4) spin model. 

The 0(4) scaling has been obtained also with improved Wilson quarks. Figure Q] 
(right) is the result for the case of the clover-improved Wilson quark action coupled 
with the RG gauge action, obtained on a 16 3 x 4 lattice. 1 ** The vertical axis is 
and the horizontal axis is i/Zi 1 /^ 5 , where j3 and 5 are the critical exponents obtained 
in the 0(4) spin model. The dashed line is the 0(4) scaling function. They fitted 
the data to the scaling function adjusting /3 ct and the scales of two axes. As seen 
from Fig. [1] (right), QCD data is well described by the 0(4) scaling ansatz. This 
result suggests that the chiral phase transition is of second order for 2 flavor QCD. 



§4. Critical temperature of 2 flavor QCD 



The critical temperature (T c ) is one of the most fundamental quantities in the 
QCD thermodynamics and is important in phenomenological studies of heavy ion 
collisions. Recently several groups 14 ^ 16 - 1 have tried to determine T c near the physical 
mass parameter in 2+1 flavor QCD by simulations with improved staggered quarks. 
However, the results are still contradictory to each other. 

The WHOT-QCD Collaboration reported recently a preliminary result of T c for 
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Fig. 2. Left: The pseudo-critical point (3 pc as a function of m q ~ — 1/A" C for JVt = 4 (circle) 
and Nt = 6 (square). 4 ' Right: Comparison of T pc scaled by ro between the staggered quark 
action (open symbol) 15 ' and the Wilson quark action (filled symbol) 4 - 1 for Nt = 4 and 6. 



2 flavor QCD renewing the analysis done in Ref. I). 4 -* They determined the pseudo- 
critical points (3p C defined from the peak of the Polyakov loop susceptibility on 16 3 x 4 
and 16 3 x 6 lattices, as a function of the hopping parameter K. 

As seen in the previous section, the subtracted chiral condensate satisfies the 
scaling behavior with the critical exponents and scaling function of the 3-dimensional 
0(4) spin model. Assuming that the pseudo-critical temperature from the Polyakov 
loop susceptibility follows the same scaling law as the 0(4) spin model, i.e. t pc ~ h y 
with y = P6 = 0.537(7), the data of (3 pc in Fig. [2] (left) are fitted by j3 pc = pa+Ah 1 ^ 
with two free parameters, ft c t and A, in the range of /3 = 1.8-1.95 for Nt = 4 and 
P = 1.95-2.10 for Nt = 6. First, they adopted the definition m q a ~ l/K - 1/K C 
as the quark mass where K c is the chiral point where the pion mass vanishes at 
T = for each /3. The critical temperature T c is calculated in the chiral limit 
using T = l/(Nfa). The lattice spacing a is estimated from the vector meson mass 
assuming my(T = 0) = m p = 770 MeV at [3 c t on K c . By this procedure, They 
obtained preliminary results of T c = 183(3) MeV for N = 4 and 174(5) MeV for 
N t = 6. They also calculated (3^ using the relation of m^ WI oc rrtpg, where mps is 
the pseudo-scalar meson mass and m^ WI is the quark mass obtained from the axial 
vector Ward-Takahashi identity. The results of T c are 173(3) MeV (N t = 4), 167(3) 
MeV (JV t = 6) for h = (m PS a) 2 and 176(3) MeV (N t = 4) for h = m^ WI a. It is 
noted that these 0(4) fits reproduce the data of /3 pc much better than a linear fit 
/3 pc = f3 c t + Ah. A tentatively conclusion is that the critical temperature in the chiral 
limit is in the range 170-186 MeV for N t = 4 and 164-179 MeV for N t = 6. There 
is still a large uncertainty from the choice of the fit ansatz. To remove this, further 
simulations at lighter quark masses are necessary. 

Next, we compare these results with those of a staggered quark action. We 
plot the results of the pseudo-critical temperature (T pc ) in unit of Sommer scale (ro) 
as a function of mps^o in Fig. [2] (right) together with those by the RBC-Bielefeld 
Collaboration using 2+1 flavor p4-improved staggered quark action. 15 ) As seen in 
this figure, results of T pc obtained by different quark actions seem to approach the 
same function of rra P sro as Nt increases. 
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§5. Equation of state at /x q = and fi q ^ 

The studies of the equation of state (EoS) can provide basic input for the analysis 
of the experimental signatures for QGP formation, e.g. the EoS will control the 
properties of any hydrodynamic expansion. 

For the study at zero chemical potential, the integral method is commonly used. 
This method is based on the equation for pressure, p = (T /V)lnZ , where Z is 
the partition function. Because the derivatives of the partition function can be 
expressed by expectation values of operators, which are computable by a Monte- 
Carlo simulation, we obtain the pressure by integrating this expectation value in the 
parameter space. For the case of the Wilson quark, we have 

t* - - n < {-m (f ) - <vaiue at t = o) i (54) 

with d£ = (df3',dK r ) on the integration path. The starting point of the integration 
path should be chosen such that p ~ there. 

The CP-PACS Collaboration carried out a systematic calculation of EOS with a 
2 flavor Wilson-type quark action at fj, q = 0. They used the RG gauge and the clover 
quark actions with 16 3 x 4 and 16 3 x 6 lattices. 2 ) The translation from the results 
obtained by Eq. (|5-ip as functions of K) to those of physical parameters can be 
done using Fig. [3] (left). The thin solid lines shows the lines of constant physics 
(LCP), which they determine by mps/my (the ratio of pseudo-scalar and vector 
meson masses at T = 0). The chiral limit K c corresponds to LCP for mpg/my = 0. 
The bold dashed line denoted as K t (N t = 4) represents the pseudo-critical line 
T/T pc = 1 at Nf = 4. The thin dashed lines represent the lines of constant T/T pc 
estimated by T jray = (A^mya) -1 . 

The right panel of Fig. [3] is the result of the pressure as a function of temperature 
for each LCP. Filled and open symbols are the results for Nt = 4 and 6 respectively. 
Different shapes of the symbol correspond to different values of mps/mv, i.e. differ- 
ent quark masses. This figure shows that the pressure is almost independent of the 
quark mass in a wide range of mps/rny. However, the A^-dependence is sizeable, 



S. Ejiri 



• in 2 



SB (N=4)- 



SB limit 



1.5 

277 



• d 2 (x q /^)/d(ji q /D 2 



/I d □ □ S □ 



SB (N =4) 

T 

SB limit 



1.5 

77T 



Fig. 4. Left: Quark number (circle) and isospin (square) susceptibilities at fi q — hi — 0. Right: 
The second derivatives of these susceptibilities. 



hence further simulations with large Nt are important. 

On the other hand, the studies of the EoS at non-zero baryon number density is 
attracting interest widely since the chemical potential dependence of the EoS should 
explain the difference of experimental results obtained at different beam energy. 
Moreover, hadronic fluctuations at finite densities are closely related to the appear- 
ance of the critical point in the (T, n q ) plane and may be experimentally testable 
by an event-by-event analysis of heavy ion collisions. The fluctuations can also be 
studied by numerical simulations of lattice QCD calculating the quark number and 
isospin susceptibilities, Xq an d Xi- They correspond to the second derivatives of the 
pressure with respect to \x q and where \xi is the isospin chemical potential. From 
a phenomenological argument in the sigma model, Xq is singular at the critical point, 
whereas xi shows no singularity there. 

To investigate the EoS at finite density, the Bielefeld-Swansea Collaboration 
proposed a Taylor expansion method. 17 ** Since Monte-Carlo simulations cannot be 
performed at finite density due to the sign problem, one evaluates the \x q dependence 
of InZ by computing the higher order derivatives of lnZ instead of integrating the 
first derivative. 

The WHOT-QCD Collaboration performed simulations at mpg/my = 0.65 and 
0.80 on a 16 3 x 4 lattice with the improved Wilson quark action. They calculated 
the second and forth derivatives of pressure which correspond to the Xq an d Xi an d 
their second derivatives with respect to [i q and \xi at fj, q = fij = 0. (Note that the 
odd derivatives are zero at fx q = 0.) 

The left panel of Fig. [Ushows Xq/T 2 (circle) and xi/T 2 (square) at mps/my = 
0.8 and fj, q = fij = as functions of T/T pc . It is found that Xq/T 2 and xi/T 2 increase 
sharply at T pc , in accordance with the expectation that the fluctuations in the QGP 
phase are much larger than those in the hadron phase. Their second derivatives 
d 2 (xq/T 2 )/d(fi q /T) 2 and d 2 (xi/T 2 )/d(fi q /T) 2 are shown in Fig. H (right). The basic 
features are quite similar to those found previously with the p4-improved staggered 
fermions. 17 ) d 2 (xi/T 2 )/d(/j q /T) 2 remains small around T pc , suggesting that there 
are no singularities in xi at non-zero density. On the other hand, we expect a large 
enhancement in the quark number fluctuations near T pc as approaching the critical 
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point in the (T, plane. The dashed line in Fig. [4] (right) is a prediction from the 
hadron resonance gas model at low temperature, d 2 Xq/d^L 2 ~ 9x q /T 2 . Although 
current statistical errors in Fig. H (right) are still large, d 2 (x q /T 2 )/d(fi q /T) 2 near 
T pc is much larger than that at high temperatures. At the right end of the figure, 
values of free quark-gluon gas (Stefan-Boltzmann gas) for Nt = 4 and for Nt = oo 
limit are shown. Since the lattice discretization error in the EoS is known to be large 
at Nt = 4 with their quark action, it is needed to extend this study to larger Nt for 
the continuum extrapolation. 

§6. Heavy quark free energies 

Finally, we discuss a free energy between static quarks. Clarification of the 
interaction between heavy quarks in QGP is important to understand the properties 
of charmoniums in heavy ion collisions. The heavy quark free energy Fm in various 
color channels M can be measured separately on the lattice by the correlations of 
the Polyakov loop with an appropriate gauge fixing. 

Recently, the heavy quark free energy for 2 flavor QCD with dynamical Wilson 
quarks in the Coulomb gauge are studied by the WHOT-QCD Collaboration, using 
the same configuration for the calculation of the EoS in the previous section. 5 ) With 
the improved actions they adopted, the rotational symmetry is well restored in the 
heavy quark free energies, 18 - 1 hence it is not necessary to introduce terms correcting 
lattice artifacts at short distances to analyze the data. 

They found that, at T > T pc , the free energies of QQ and QQ normalized to 
be zero at large separation show attraction (repulsion) in the color singlet and anti- 
triplet channels (color octet and sextet channels), and fitted the free energy data in 
each channel by the screened Coulomb form, 

Fm (?, T) — Fm (oo, T) = C(M)^^- e -™^ T > , (6-1) 

r 

where a e g(T) and mo{T) are the effective running coupling and Debye screening 
mass, respectively. The Casimir factor C(M) = (X^a=i ' for color channel 
M is explicitly given by C(l) = -|,C(8) = g,C(6) = §,C(3*) = -§. The results 
of a e fr(T) and mz>(T) are shown in Fig. [5] for mps/my = 0.65. They found that 
there is no significant channel dependence in a e g(T) and mx>(T) at sufficiently high 
temperatures (T <; 2T pc ). In other words, the channel dependence in the free energy 
can be well absorbed in the kinematical Casimir factor at high temperatures, 

The magnitude and the T-dependence of mr){T) is consistent with the next- 
to-leading order calculation in thermal perturbation theory. Moreover it is also 
well approximated by the leading order form with an "effective" running coupling 
defined from a e s{T). On the other hand, by comparing these results with the results 
by an improved staggered quark action, 19 ) they found that a e g(T) does not show 
appreciable difference while mi)(T) in the Wilson quark action is larger than that 
of the staggered quark action by 20%. To draw a definite conclusion, however, 
simulations with smaller lattice spacings, i.e., larger lattice sizes in the temporal 
direction (such as Nt = 6 or larger) at smaller quark masses are required. 
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Fig. 5. The effective running coupling a e ff(T) (left) and Debye screening mass m_o(T) (right) for 
each color channel as a function of temperature at mps/my — 0.65. 



§7. Sumarry 

We reported the current status of the study of QCD thermodynamics using a 
Wilson quark action. The basic properties of QCD with 2 flavors of Wilson quarks 
such as the phase structure and the scaling behavior have been already studied. 
Moreover, the qualitative agreements between the results obtained by a Wilson quark 
action and a staggered quark action have been discussed for some thermodynamic 
quantities. The quantitative improvements are important in the future studies. 

I would like to thank the members of the WHOT-QCD Collaboration for col- 
laboration, discussions and comments. This work is supported by Grants-in-Aid of 
the Japanese MEXT (No. 18740134) and Sumitomo Foundation (No. 050408). 
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